#include "plugin.h"

#include <assert.h>
#include <GL\glew.h>

#define min(a,b)            (((a) < (b)) ? (a) : (b))
#define max(a,b)            (((a) > (b)) ? (a) : (b))

// --- forward declarations ---------------------------------------------------

void ProcessItem(BWInstanceData* instance, LWItemID itemID);
void ProcessPolygon(BWInstanceData* instance, BWcoInstanceData* coInstance, 
                    const PolygonInfo* polInfo);
void ProcessTriangle(BWInstanceData* instance, BWcoInstanceData* coInstance, 
                     const PolygonInfo* triangle);
void PrintPolygonInfo(const PolygonInfo* polInfo);
int CreateParticlesInTriangle(
    BWInstanceData* instance, BWcoInstanceData* coInstance, 
    const PolygonInfo* triangle);
void RenderTriangleInZBuffer(
    const double wp1[3], const double wp2[3], const double wp3[3]);
size_t PolScan(void* data, LWPolID polID);
BOOL IsTriangleVisible(const double p1[3], const double p2[3], 
                       const double p3[3], double frameAspect, double margin);
BOOL IsVisible(const double screenCoords[3], double frameAspect, double margin);
double GetFrameAspect(const LWSceneInfo* sceneInfo);
void GetViewVector(const BWInstanceData* instance, const double* point, 
                   double* viewVector);
void DrawTriangle(const double wp1[3], const double wp2[3], const double wp3[3],
                  float r, float g, float b);

typedef struct st_PolScanData
{
    BWInstanceData* bwInstance;
    BWcoInstanceData* bwCoInstance;
    LWMeshInfo* mesh;    
} PolScanData;

// ----------------------------------------------------------------------------
void ProcessScene(BWInstanceData* instance)
{   
    int totalNumPoints;
    int countPoints = 0;
    LWItemID itemID, childItemID;
    const LWItemInfo* itemInfo = getLWItemInfo(instance->globalFunc);
    const LWObjectInfo* objectInfo = getLWObjectInfo(instance->globalFunc);

    // triangle count will act as random seed
    instance->numProcessedTriangles = 0;
    instance->projMatrixDirty = TRUE;
   
    // make sure there are objects in the scene
    totalNumPoints = getLWSceneInfo(instance->globalFunc)->numPoints;
    if (totalNumPoints <= 0) {
        WriteTrace("No objects in the scene, doing nothing");
        return;
    } else {
        WriteTrace("Found %d points", totalNumPoints);
    }

    // OpenGL drawing is used as a scratch buffer to find polygon edges
    SetupGLMatrices(instance->globalFunc, instance->filterAccess->width, 
                    instance->filterAccess->height);

    // adjust OpenGL settings
    glShadeModel(GL_FLAT);							    
    glClearColor(0.0f, 0.0f, 0.0f, 0.0f); // Black Background    
    glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST); 
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    
    glDisable(GL_LIGHTING);
    glDisable(GL_BLEND);
    glDisable(GL_TEXTURE_2D); 
    glDisable(GL_DEPTH_TEST);

    itemID = itemInfo->first(LWI_OBJECT, LWITEM_NULL);
    while (itemID != LWITEM_NULL) {
        unsigned int itemFlags;
        
        WriteTrace("Found object '%s' (itemID: %d), containing %d points.", 
            itemInfo->name(itemID), itemID, objectInfo->numPoints(itemID));
        
        itemFlags = itemInfo->flags(itemID);
        if (!(itemFlags & LWITEMF_ACTIVE))
        {
            WriteTrace("  Ignoring object (it's inactive).");
        }
        else
        {        
            if (objectInfo->numPoints(itemID) >= 3)
            {   
                countPoints += objectInfo->numPoints(itemID);
                ProcessItem(instance, itemID);        
                childItemID = itemInfo->firstChild(itemID);        
                if (childItemID != LWITEM_NULL)
                    WriteTrace("  Found child item: %d", childItemID);
            }
        }

        itemID = itemInfo->next(itemID);
    }

    checkGLerror();
}

// ----------------------------------------------------------------------------
void SetupGLMatrices(GlobalFunc* globalFunc, int width, int height)
{
    LWItemID camID;
    LWTime time;
    double hfov, vfov, zoomFactor;
    double r[3], u[3], f[3], p[3];

    const LWItemInfo* itemInfo = getLWItemInfo(globalFunc);
    const LWSceneInfo* sceneInfo = getLWSceneInfo(globalFunc);
    const LWCameraInfo* camInfo = getLWCameraInfo(globalFunc);
    const LWTimeInfo* timeInfo = getLWTimeInfo(globalFunc);

    // --- define viewport ---
    glViewport(0, 0, width, height);

	// --- adjust projection matrix ---
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();	
	
    time = timeInfo->time;
    camID = sceneInfo->renderCamera(time);

    camInfo->fovAngles(camID, time, &hfov, &vfov);
    zoomFactor = camInfo->zoomFactor(camID, time);

    gluPerspective(rad2deg(vfov), GetFrameAspect(sceneInfo), 0.1, 1000);
    
    // --- adjust modelview matrix ---

    itemInfo->param(camID, LWIP_RIGHT,      time, r);
    itemInfo->param(camID, LWIP_UP,         time, u);
    itemInfo->param(camID, LWIP_FORWARD,    time, f);
    itemInfo->param(camID, LWIP_W_POSITION, time, p);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    // up vector is inverted, since openGL uses a right-handed coordinate-
    // system, while Lightwave's is left-handed
    gluLookAt(
        p[0], p[1], p[2], 
        p[0] + f[0], p[1] + f[1], p[2] + f[2],
        -u[0], -u[1], -u[2]);
}

// ----------------------------------------------------------------------------
void ProcessItem(BWInstanceData* instance, LWItemID itemID)
{
    BWcoInstanceData* coInstance;
    LWMeshInfo *mesh;
    PolScanData scanData;
    
    coInstance = GetFromGlobalList(itemID);
    if (!coInstance)
    {
        WriteTrace("Object %d contains no brushwave settings -> ordinary rendering", itemID);
    }
    else
    {
        WriteTrace("Using the following settings for object %d:", itemID);
        WriteTrace("  Brush Texture: %s", coInstance->renderSettings.brushFilename);
        WriteTrace("  Particle Size: %f", coInstance->renderSettings.partSize);
        WriteTrace("  Particle Density: %f", coInstance->renderSettings.partDensity);
        WriteTrace("  Colour Variance: %f", coInstance->renderSettings.colorVariation); 
        WriteTrace("  Colour Displacement: %f", coInstance->renderSettings.colorDisplacement);
        WriteTrace("  Align Particles: %d", coInstance->renderSettings.alignParticles);
        WriteTrace("  Backface Culling: %d", coInstance->renderSettings.backFaceCulling);
    }

    mesh = getLWObjectInfo(instance->globalFunc)->meshInfo(itemID, 1);

    scanData.bwInstance = instance;
    scanData.bwCoInstance = coInstance;
    scanData.mesh = mesh;    
    
    mesh->scanPolys(mesh, PolScan, &scanData);
}

// ----------------------------------------------------------------------------
size_t PolScan(void* data, LWPolID polID)
{    
    int i;
    LWID polType;
    PolygonInfo polInfo;    
    PolScanData* scanData = (PolScanData*) data;
    BWInstanceData* instance = scanData->bwInstance;
    BWcoInstanceData* coInstance = scanData->bwCoInstance;
    LWMeshInfo* mesh = scanData->mesh;
    LWFVector normal;

    polInfo.numPoints = mesh->polSize(mesh, polID);
    polType = mesh->polType(mesh, polID);
    
    // skip polygons with a zero-area
    if (polInfo.numPoints < 3) return 0;

    //WriteTrace("Polygon %d, type %d -> numPoints: %d", polID, polType, polInfo.numPoints);

    polInfo.normal = calloc(3, sizeof(double));
    mesh->polOtherNormal(mesh, polID, normal);

    polInfo.normal[0] = (double) normal[0];
    polInfo.normal[1] = (double) normal[1];
    polInfo.normal[2] = (double) normal[2];

    polInfo.points = calloc(polInfo.numPoints, sizeof(double*));

    for (i=0; i<polInfo.numPoints; ++i)
    {          
        LWFVector points;
        LWPntID pointID;        
        
        // get point and save it in polInfo
        pointID = mesh->polVertex(mesh, polID, i);
        mesh->pntOtherPos(mesh, pointID, points);

        polInfo.points[i] = calloc(3, sizeof(double));
        polInfo.points[i][0] = (double) points[0];
        polInfo.points[i][1] = (double) points[1];
        polInfo.points[i][2] = (double) points[2];
    }

    ProcessPolygon(instance, coInstance, &polInfo);

    // free created memory
    for (i=0; i<polInfo.numPoints; ++i)
        free(polInfo.points[i]);

    free(polInfo.points);
    free(polInfo.normal);

    return 0;
}

// ----------------------------------------------------------------------------
void ProcessPolygon(BWInstanceData* instance, BWcoInstanceData* coInstance,
                    const PolygonInfo* polInfo)
{
    int i, numPoints;

    numPoints = polInfo->numPoints;
    assert(numPoints >= 3);

    // polygon could have any number of vertices -> we split it up in 
    // as many triangles as needed.
    for (i=0; i<numPoints-2; ++i)
    {
       PolygonInfo triangle;
       triangle.numPoints = 3;
       triangle.points = calloc(3, sizeof(double*));

       triangle.points[0] = polInfo->points[0];
       triangle.points[1] = polInfo->points[i+1];
       triangle.points[2] = polInfo->points[i+2];

       triangle.normal = polInfo->normal;

       ProcessTriangle(instance, coInstance, &triangle);

       free(triangle.points);
    }
}

// ----------------------------------------------------------------------------
void ProcessTriangle(BWInstanceData* instance, BWcoInstanceData* coInstance,
                     const PolygonInfo* triangle)
{
    assert(triangle->numPoints == 3);
    
    if (coInstance)
    {
        // before creating particles, check if triangle can be culled!
        BOOL cullFace = FALSE;
        
        if (coInstance->renderSettings.backFaceCulling)
        {
            double viewVector[3];
            GetViewVector(instance, triangle->points[0], viewVector);
            cullFace = (dot(triangle->normal, viewVector) < 0);
        }
        
        if (!cullFace)
        {
            DrawTriangle(triangle->points[0], triangle->points[1], triangle->points[2], 
                1.0f, 1.0f, 1.0f);
            
            CreateParticlesInTriangle(instance, coInstance, triangle); 

            DrawTriangle(
                triangle->points[0], triangle->points[1], triangle->points[2], 
                0.0f, 0.0f, 0.0f);
        }
    }
    else
        RenderTriangleInZBuffer(
            triangle->points[0], triangle->points[1], triangle->points[2]); 

    instance->numProcessedTriangles++;
}

// ----------------------------------------------------------------------------
// --- this is where most of the magic happens ;) -----------------------------
// --- ------------------------------------------------------------------------
int CreateParticlesInTriangle(
    BWInstanceData* instance, BWcoInstanceData* coInstance, 
    const PolygonInfo* triangle)
{
    static int lastIX = -1, lastIY = -1;
    GLfloat pixelRegion[9];
    int i, ix, iy, colorRefPos, depthRefPos, numParticles;
    float p, cos, alpha, caNumParticles, fadingZoneFactor, 
          fadingZoneStart, fadingZoneEnd, rndAddition, colVariance,
          finalPoint[3];
    double *wp1, *wp2, *wp3, *normal,
           p1[3], p2[3], p3[3],
           tmp1[3], tmp2[3], p1p2[3], p1p3[3], 
           testVector[3], testPoint[3], wTestPoint[3],
           rnd1, rnd2, area, s, x, y, 
           minScrX, maxScrX, minScrY, maxScrY, minScrZ, maxScrZ,
           triangleScrWidth, triangleScrHeight, triangleScrDepth, zTestOffset,
           pixelScrWidth, pixelScrHeight, frameAspect;
    BOOL isInside;
    Particle* newParticle;
    RenderSettings renderSettings = coInstance->renderSettings;
    const LWSceneInfo* sceneInfo = getLWSceneInfo(instance->globalFunc);
    TextureInfo* textureInfo;

    wp1 = triangle->points[0]; // points in world coords
    wp2 = triangle->points[1];
    wp3 = triangle->points[2];
    normal = triangle->normal;

    // options:
    fadingZoneFactor = 0.2f;
    numParticles = 0;
    frameAspect = GetFrameAspect(sceneInfo);
    
    ConvertToScreenCoords(instance, wp1, p1);
    ConvertToScreenCoords(instance, wp2, p2);
    ConvertToScreenCoords(instance, wp3, p3);

    //WriteTrace("  Point1 sc: (%f, %f, %f)", p1[0], p1[1], p1[2]);
    //WriteTrace("  Point2 sc: (%f, %f, %f)", p2[0], p2[1], p2[2]);
    //WriteTrace("  Point3 sc: (%f, %f, %f)", p3[0], p3[1], p3[2]);

    if (!IsTriangleVisible(p1, p2, p3, frameAspect, 0.0))
        return 0;
    
    init_genrand(instance->numProcessedTriangles);

    // calculate some properties of the triangle we will need later
    minScrX = min(min(p1[0], p2[0]), p3[0]);
    maxScrX = max(max(p1[0], p2[0]), p3[0]);
    minScrY = min(min(p1[1], p2[1]), p3[1]);
    maxScrY = max(max(p1[1], p2[1]), p3[1]);
    minScrZ = min(min(p1[2], p2[2]), p3[2]);
    maxScrZ = max(max(p1[2], p2[2]), p3[2]);
    triangleScrWidth  = maxScrX - minScrX;
    triangleScrHeight = maxScrY - minScrY;
    triangleScrDepth  = maxScrZ - minScrZ; 
    pixelScrWidth  = 2.0 * frameAspect / sceneInfo->frameWidth;
    pixelScrHeight = 2.0 / sceneInfo->frameHeight;
    
    // zTestOffset = maximum depth difference that can occur in one
    //               pixel of the triangle    
    zTestOffset = triangleScrDepth / 
        min(triangleScrWidth / pixelScrWidth, triangleScrHeight / pixelScrHeight) * 
        2.0;
    zTestOffset = max(min(pixelScrWidth, pixelScrHeight), zTestOffset);

    /*

    // following: just a test

    // we will draw the polygon onto the pbuffer to have a reference about
    // where the actual pixels are in the image. However, the aliasing results
    // in Lightwave differ from OpenGL, so the triangle we draw in OpenGL
    // is resized a tiny little bit before drawing it.

    if (FALSE)
    {
        double center[3], c_wp1[3], c_wp2[3], c_wp3[3],
            c_wp1_scaled[3], c_wp2_scaled[3], c_wp3_scaled[3];
        double scaleVal = 0.6f;
        
        center[0] = (wp1[0] + wp2[0] + wp3[0]) / 3.0;
        center[1] = (wp1[1] + wp2[1] + wp3[1]) / 3.0;
        center[2] = (wp1[2] + wp2[2] + wp3[2]) / 3.0;

        subtract(c_wp1, wp1, center);
        subtract(c_wp2, wp2, center);
        subtract(c_wp3, wp3, center);

        multiply(c_wp1_scaled, c_wp1, scaleVal);
        multiply(c_wp2_scaled, c_wp2, scaleVal);
        multiply(c_wp3_scaled, c_wp3, scaleVal);

        add(resizedWp1, center, c_wp1_scaled);
        add(resizedWp2, center, c_wp2_scaled);
        add(resizedWp3, center, c_wp3_scaled);
    }
    
    */
    
    // create vectors that spawn up the triangle
    subtract(p1p2, wp2, wp1);
    subtract(p1p3, wp3, wp1);
    
    area = calcTriangleArea2D(p1, p2, p3);    
    // test: fixed particle count
    //area = calcTriangleArea3D(wp1, wp2, wp3);
    
    // ca = circa; 4 * frameAspect = the screen area 
    caNumParticles = (float) (area * (renderSettings.partDensity / renderSettings.partSize)); 

    //WriteTrace("Area:               %f", area);
    //WriteTrace("Expected particles: %f", caNumParticles);

    alpha = 1.0;

    // --- create random particles ---
    p = 0;
    while (alpha > 0)
    {
        rndAddition = (float) genrand_real1();
        colVariance = (float) genrand_real1() * (float) renderSettings.colorVariation * 2.0f - 
            renderSettings.colorVariation;
        
        isInside = FALSE;
        while(!isInside)
        {
            // create random point in triangle
            rnd1 = genrand_real1();
            rnd2 = genrand_real1();

            multiply(tmp1, p1p2, rnd1);
            multiply(tmp2, p1p3, rnd2);

            add(testVector, tmp1, tmp2);
            add(wTestPoint, wp1, testVector);

            isInside = isPointInTriangle(wTestPoint, wp1, wp2, wp3);
        }

        // soft fading of particles
        fadingZoneEnd = p + rndAddition;
        fadingZoneStart = (1.0f - fadingZoneFactor) * fadingZoneEnd;
        if (caNumParticles >= fadingZoneEnd) alpha = 1.0f;
        else if (caNumParticles <= fadingZoneStart) alpha = 0.0f;
        else alpha = (caNumParticles - fadingZoneStart) / 
                     (fadingZoneEnd * fadingZoneFactor);
        
        // at first glance, it would seem to be a good idea to ignore the particle
        // if it has an alpha value of zero -- however, it leads to more flickering
        // when there are many small polygons. I think the reason is that fewer 
        // "lastIX" and "lastIY" are calculated. Hence, I removed this optimization.
        //
        // if (alpha == 0) return numParticles;

        if (renderSettings.alignParticles || renderSettings.backFaceCulling)
        {
            // when particle is aligned, its screen area is approx. the
            // dot product of the view vector and the normal vector
            
            double viewVector[3];
            GetViewVector(instance, wTestPoint, viewVector);
            normalize(viewVector, viewVector);
            cos = (float) dot(normal, viewVector);
        }
        else
        {
            cos = 1.0f;
        }

        if (renderSettings.backFaceCulling && cos < 0) 
            return 0;
        
        if (!renderSettings.alignParticles) 
            cos = (cos < 0 ? -1.0f : 1.0f);

        p += fabsf(cos);

        // particle position in screen coordinates
        ConvertToScreenCoords(instance, wTestPoint, testPoint);
        
        // ... and in actual pixels
        s = sceneInfo->frameHeight * 0.5;
        x = sceneInfo->frameWidth * 0.5 + testPoint[0] * s / sceneInfo->pixelAspect;
        y = s - testPoint[1] * s;

        ix = (int) x;
        iy = (int) y;

        // check new xy-bounds (must be within image size)
        if (ix < 0 || ix >= instance->filterAccess->width ||
            iy < 0 || iy >= instance->filterAccess->height) continue;
        
        // reduce flickering: check if particle position is inside 
        // aliased drawing of triangle
        ReadPixels(ix-1, iy-1, 3, 3, GL_RED, pixelRegion);

        if (pixelRegion[4] == 0) 
        {
            // center of particle is rendered OUTSIDE of triangle
            // -> this could lead to flickering, so we change the position
            
            int i, j, found_i, found_j;
            found_i = found_j = -1;
            
            for (i=0; i<3; ++i)
            {
                for (j=0; j<3; ++j)
                {
                    float pixVal = pixelRegion[3*i+j];
                    if (pixVal > 0.9f)
                    {
                        // pixel in current triangle
                        found_i = i;
                        found_j = j;
                        break;
                    }
                }
                if (found_i != -1) break;
            }

            if (found_i != -1)
            {
                // found pixel in surrounding area
                lastIX = ix = ix - 1 + found_j;
                lastIY = iy = iy - 1 + found_i;
            }
            else
            {
                // no suitable pixel found :(
                if (lastIX == -1)
                    continue;
                else
                {
                    // use the last working particle position.
                    // that looks better than kicking the particle completely.
                    ix = lastIX;
                    iy = lastIY;
                }
            }
            
            // check new xy-bounds (must be within given reference picture)
            if (ix < 0 || ix >= instance->filterAccess->width ||
                iy < 0 || iy >= instance->filterAccess->height) continue;
        }
        else
        {
            lastIX = ix;
            lastIY = iy;
        }

        // prepare texture
        textureInfo = GetTextureInfo(instance, coInstance->renderSettings.brushFilename);
        if (!textureInfo) return 0; // something went wrong! panic and return ;)
        
        // if we are here, we found a valid particle - so let's create it.
        for (i=0; i<3; ++i) finalPoint[i] = (float) wTestPoint[i];
        
        newParticle = CreateParticle(finalPoint, testPoint[2], 
            renderSettings.partSize, renderSettings.alignParticles);
        
        newParticle->textureInfo = textureInfo;
        newParticle->aspectRatio = (float) textureInfo->width / 
            (float) textureInfo->height;

        // get position for access of reference images
        depthRefPos = (instance->filterAccess->width * iy + ix);
        colorRefPos = depthRefPos * 4;        
        
        // save normal vector
        for (i=0; i<3; ++i) newParticle->normal[i] = (float) normal[i];

        if (testPoint[2] <= instance->depthBuffer[depthRefPos] + zTestOffset)
        {        
            float newColor[4], greyValue, offset;
            for (i=0; i<3; ++i) newColor[i] = instance->pixelBuffer[colorRefPos + i] + colVariance;
            newColor[3] = alpha;
            
            // Z Color Displacement
            if (renderSettings.colorDisplacement != 0)
            {
                greyValue = 0.3f * newColor[0] + 0.59f * newColor[1] + 0.11f * newColor[2];
                offset = (float) (renderSettings.partSize * newParticle->zDist * (1.0 - greyValue));
                newParticle->colorDisplacement = renderSettings.colorDisplacement * offset;
            }
            else
            {
                newParticle->colorDisplacement = 0;
            }

            // colour of particle is received from first renderstep (Lightwave)
            for (i=0; i<3; ++i) newParticle->colour[i] = newColor[i];
            newParticle->colour[3] = alpha;
            
            numParticles++;
            instance->particleList = InsertFront(instance->particleList, newParticle);
        }
        else
        {            
            // display particle in red (testing purposes)
            newParticle->colour[0] = 1;
            newParticle->colour[1] = 0;
            newParticle->colour[2] = 0;
            newParticle->colour[3] = 1;

            //instance->particleList = InsertFront(instance->particleList, newParticle);
            free(newParticle);
            
            // ignore particle, since it is on a probably invisible triangle!
             //WriteTrace("Invisible particle! z = %f, zDiff = %f", 
             //   newParticle->pos[2],
             //   newParticle->pos[2] - depthBuffer[depthRefPos]);
        }
    }

    return numParticles;
}

void RenderTriangleInZBuffer(
    const double wp1[3], const double wp2[3], const double wp3[3])
{        
    glPushAttrib(GL_ALL_ATTRIB_BITS);
    
    // only draw in depth buffer
    glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE); 
    glDepthMask(GL_TRUE);
    glEnable(GL_DEPTH_TEST);

    glBegin(GL_TRIANGLES);		
        glVertex3dv(wp1);
        glVertex3dv(wp2);    
        glVertex3dv(wp3);
    glEnd();    

    glPopAttrib();
}

// ----------------------------------------------------------------------------
// returns coordinates within { [-frameAspect, frameAspect], [-1, 1] }
// ----------------------------------------------------------------------------
void ConvertToScreenCoords(BWInstanceData* instance, 
                           const double* worldCoords, double* screenCoords)
{
    LWItemID camID;
    LWTime time;
    static double u[3], v[3], n[3], r[3], e, w;
    static MAT4 t;
    double zdist;
    const LWItemInfo* itemInfo = getLWItemInfo(instance->globalFunc);
    const LWSceneInfo* sceneInfo = getLWSceneInfo(instance->globalFunc);
    const LWCameraInfo* camInfo = getLWCameraInfo(instance->globalFunc);
    const LWTimeInfo* timeInfo = getLWTimeInfo(instance->globalFunc);
    
    time = timeInfo->time;
    camID = sceneInfo->renderCamera(time);
    itemInfo->param(camID, LWIP_RIGHT,        time, u);
    itemInfo->param(camID, LWIP_UP,           time, v);
    itemInfo->param(camID, LWIP_FORWARD,      time, n);
    itemInfo->param(camID, LWIP_W_POSITION,   time, r);
    e = -1.0 * camInfo->zoomFactor(camID, time);

    if (instance->projMatrixDirty)
    {
        instance->projMatrixDirty = FALSE;

        t[0][0] = u[0];
        t[0][1] = v[0];
        t[0][2] = n[0];
        t[1][0] = u[1];
        t[1][1] = v[1];
        t[1][2] = n[1];
        t[2][0] = u[2];
        t[2][1] = v[2];
        t[2][2] = n[2];
        t[3][0] = -dot( r, u );
        t[3][1] = -dot( r, v );
        t[3][2] = -dot( r, n ) + e;
        t[0][3] = -n[0] / e;
        t[1][3] = -n[1] / e;
        t[2][3] = -n[2] / e;
        t[3][3] = 1.0 - t[3][2] / e;
    }

    w = transform(worldCoords, t, screenCoords);
    
    // calculate real z distance
    zdist = w * screenCoords[2] - e;
    screenCoords[2] = zdist;
}

// ----------------------------------------------------------------------------
void ConvertMatrix(const MAT4 in, double out[16])
{
    out[0] =  in[0][0];
    out[1] =  in[1][0];
    out[2] =  in[2][0];
    out[3] =  in[3][0];
    out[4] =  in[0][1];
    out[5] =  in[1][1];
    out[6] =  in[2][1];
    out[7] =  in[3][1];
    out[8] =  in[0][2];
    out[9] =  in[1][2];
    out[10] = in[2][2];
    out[11] = in[3][2];
    out[12] = in[0][3];
    out[13] = in[1][3];
    out[14] = in[2][3];
    out[15] = in[3][3];
}

// ----------------------------------------------------------------------------
BOOL IsTriangleVisible(const double p1[3], const double p2[3], 
                       const double p3[3], double frameAspect, double margin)
{
    if (IsVisible(p1, frameAspect, margin) &&
        IsVisible(p2, frameAspect, margin) &&
        IsVisible(p3, frameAspect, margin))
    {
        return TRUE;
    }

    // check regions ...
    
    if (   // above screen
        ((p1[1] - margin > 1) && (p2[1] - margin > 1) && (p3[1] - margin > 1))
        || // below screen
        ((p1[1] + margin < -1) && (p2[1] + margin < -1) && (p3[1] + margin < -1))
        || // left of screen
        ((p1[0] + margin < -frameAspect) && (p2[0] + margin < -frameAspect) && (p3[0] + margin < -frameAspect))
        || // right of screen
        ((p1[0] - margin > frameAspect) && (p2[0] - margin > frameAspect) && (p3[0] - margin > frameAspect))
        )
    {
        return FALSE;
    }

    // we are not really sure yet; exact tests could, however, be more 
    // expensive than just calculating the particles within this polygon.
    return TRUE;
}

// ----------------------------------------------------------------------------
BOOL IsVisible(const double screenCoords[3], double frameAspect, double margin)
{
    // the point is behind the camera ...
    if (screenCoords[2] <= 0) return FALSE;

    // Y-coordinate is outside image
    if (fabs(screenCoords[1]) - margin > 1) return FALSE;

    // check X-coordinate
    if (fabs(screenCoords[0]) - margin > frameAspect) return FALSE;
      
    return TRUE;
}

// ----------------------------------------------------------------------------
double GetFrameAspect(const LWSceneInfo* sceneInfo)
{
    return sceneInfo->pixelAspect * sceneInfo->frameWidth / 
        sceneInfo->frameHeight;
}

// ----------------------------------------------------------------------------
void GetViewVector(const BWInstanceData* instance, const double* point, 
                   double* viewVector)
{
    LWItemID camID;
    LWTime time;            
    double camPos[3];
    const LWItemInfo* itemInfo = getLWItemInfo(instance->globalFunc);
    const LWTimeInfo* timeInfo = getLWTimeInfo(instance->globalFunc);
    const LWSceneInfo* sceneInfo = getLWSceneInfo(instance->globalFunc);

    time = timeInfo->time;
    camID = sceneInfo->renderCamera(time);
    itemInfo->param(camID, LWIP_POSITION, time, camPos);

    subtract(viewVector, camPos, point);        
}

// ----------------------------------------------------------------------------
void DrawTriangle(const double wp1[3], const double wp2[3], const double wp3[3],
                  float r, float g, float b)
{
    glColor3f(r, g, b);
    glBegin (GL_TRIANGLES);		
        glVertex3dv(wp1);
        glVertex3dv(wp2);    
        glVertex3dv(wp3);
    glEnd();
}

// ----------------------------------------------------------------------------
void PrintPolygonInfo(const PolygonInfo* polInfo)
{
    int i;

    WriteTrace("Number of points: %d", polInfo->numPoints);
    for (i=0; i<polInfo->numPoints; ++i) {
        WriteTrace(" - %d: (%f, %f, %f)", i, 
            polInfo->points[i][0], polInfo->points[i][1], polInfo->points[i][2]);
    }
}